Graphene, which is a single layer of carbon atoms arranged in the honeycomb lattice, has been intensely studied recently.^[1–4] While most interest has focused on the unusual properties of Dirac fermions, the honeycomb geometry also results in an exotic state—i.e., the chiral d-wave superconductor.^[5,6] The chiral d-wave state is a topological superconducting state. Its topological invariant is a Chern number with the value C = 2. Under the open condition, Majorana edge states that traverse the gap appear, which may have potential applications in topological computations. However, this state has never been observed experimentally, even though its possible realization in graphene has attracted great interest.

A large number of theoretical works have studied the pairing mechanism in doped graphene. In the mean-field theory for the t–J model, a d-wave solution is found to be dominant over a large range of doping.^[7–9] A variational Monte Carlo study of the Hubbard Hamiltonian with intermediate interaction yields a similar result.^[10] The prediction for the low doping region is further supported by a variational study using Grassman tensor product state on the t–J model^[11] and a determinant quantum Monte Carlo (DQMC) study on the Hubbard model.^[12] Nevertheless, a functional renormalization calculation on a lightly doped honeycomb lattice finds that there is no d + id pairing instability in the simple Hubbard model and a finite antiferromagnetic Heisenberg interaction is needed for its realization.^[13] The dominant pairing symmetry by the variational cluster approximation in the vicinity of half filling is a triplet p-wave symmetry or Kekule superconductivity.^[14,15] A random phase approximation method shows that a singlet d + id is the leading superconducting pairing symmetry, while an additional staggered potential changes the favorable paring to the next-nearest-neighbor (NNN) triplet f-wave.^[16]

The band structure of the honeycomb lattice has a van Hove singularity (VHS) at the filling ρ=3/4. The interaction effect is greatly enhanced near the VHS due to the divergent density of states, which means that superconductivity may probably happen. A renormalization group study on the Hubbard model finds that the chiral d-wave state dominates over all competing orders at the VHS.^[17] Including a finite NNN hopping, a functional renormalization group method finds that a NNN d + id pairing instability is dominant at the VHS for the realistic and then the spin density wave (SDW) becomes dominant for . The critical interaction for the SDW phase decreases when the NNN hopping is reduced, and the SDW already wins for when the NNN hopping vanishes.^[18] Away from the VHS, the system still favors the d + id state. When a large enough NNN repulsion is included, a second-nearest-neighbor triplet f-wave pairing becomes competitive.^[18] However the singular-mode renormalization group and variational Monte Carlo calculations find that the system at the VHS is a chiral SDW for U=3.6t (t=2.8 eV). Then, a nearest-neighbor (NN) d + id paring becomes the leading instability when the system is doped away from the VHS, which is stable even in the presence of the NNN interaction.^[19] This result of the VHS agrees with calculations using a combination of exact diagonalization, density matrix renormalization group, variational Monte Carlo method, and quantum field theories on the Hubbard model.^[20] In particular, a recent DQMC study finds that the singlet d + id wave is the dominant pairing channel both at and away from the VHS filling for U=2t. It also observes that the chiral SDW is similarly relevant at the VHS, but weakens quickly upon doping away from the filling.^[21] However, since the two instabilities are not calculated on an equal footing, the leading instability is undetermined. The dynamical cluster approximation (DCA) calculations at the VHS show that while the chiral d + id pairing dominates at weak coupling, there is an enhanced tendency towards p-wave upon increasing the interactions.^[22]

Several studies have considered the low-filling levels. For example, a DQMC study on the Hubbard model at the filling ρ=0.4 finds that the d + id pairing is dominating for weak coupling U=2t.^[21] By including additional NNN hopping to flatten the band, a triplet chiral p-wave state is reported at ρ=0.2 for U=3t.^[23]

It is clear that the studies concerning the leading order and pairing symmetry in a doped honeycomb lattice are not consistent. This is due to the absence of precise methods for 2D interacting Hamiltonians, while each of the above-used methods has shortcomings.^[24–28] The DQMC method is a numerically-exact and unbiased approach to simulate quantum many-body models on large-scale lattices. Although the fermion-sign problem usually exists in doped systems, it can be circumvented to some extent by averaging the calculating quantity over the sign at the expense of much longer runs.^[29] There are mainly two studies applying the DQMC method to the Hubbard model describing the doped graphene,^[12,21] where results of specific fillings and weak couplings are presented. In this paper, we extend the previous studies and perform a systematic DQMC study of the pairing symmetry in a doped honeycomb lattice. We consider a full range of fillings from ρ=0 to 1 for both weak and strong couplings. Sizable lattices up to 512 sites are used, and the temperature is lowered to the extent beyond which the DQMC method fails to generate reliable results. Our systematic DQMC results provide new insights into the problem of superconductivity in A doped honeycomb lattice.

We consider a Hubbard Hamiltonian describing interacting Dirac fermions on a honeycomb lattice

The honeycomb lattice has a two-site unit cell. By a Fourier-transformation, the energy spectrum without interaction is directly obtained as with ( ). This noninteracting system is a semi-metal with two inequivalent Dirac points at . When the system is doped to 3/4 filling, the Fermi surface forms a perfect hexagon with the vortexes located at the saddle points (M points of the Brillouin zone). The density of states is logarithmically divergent at the filling, which is known as the VHS.

The Hubbard model (1) is solved numerically by means of the DQMC method.^[30–32] In DQMC, one decouples the on-site interaction term through the introduction of an auxiliary Hubbard–Stratonovich field, which is integrated out stochastically. The only errors are those associated with the statistical sampling, the finite spatial lattice size, and the inverse temperature discretization. These errors are well controlled in the sense that they can be systematically reduced as needed, and further eliminated by appropriate extrapolations. Since the “sign” problem is free at half-filling, the simulations can be carried out at very low temperatures. In addition to half-filling, the infamous sign problem^[33–35] can become severe upon lowering the temperature and increasing the interaction strength.^[36] To obtain reliable results, we control the average sign to be better than 0.5 by choosing proper temperatures. In the following simulations, we use the inverse temperature discretization , which has been verified to produce stable results. The temperature accessed in the simulations is down to T=1/12. The lattice has N=2×L ×L sites with L up to 16.

To determine the dominating pairing symmetry, the quantity we calculate is the uniform pairing susceptibility, which is defined as^[37,38]

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. The honeycomb lattice and phases of the bonds for the considered pairing symmetries. For a triplet pairing, there is an additional sign when the pairing is along the opposite direction of the arrow.

The superconducting pairing induced by the repulsive interaction has to be nonlocal, and here we consider the NN ones. The crystal symmetry group of the honeycomb lattice is with . Its irreducible representations classify the paring states. We can have s-wave, -wave, d_xy-wave singlet states and p_x-wave, p_y-wave, f-wave triplet states (see Fig. 1). Since and d_xy are the two basis functions of the representation , they are degenerate. Similarly, the two p-wave states are also degenerate (they belong to the representation ). From our finite lattice DQMC results, where any linear combination of the degenerate states has the same effective susceptibility, we can infer only that a chiral d + id (p+ip) symmetry is a candidate phase. A qualitative argument in favor of the chiral phase is that it allows a non-trivial solution of the gap equation while leaving the gap everywhere large. This suggests that it is energetically favored.

Figure 2(a) shows the effective susceptibility as a function of ρ at U=2t. While the values of the p- and f-channels are negative, that of the d-wave symmetry is always positive. This implies that the corresponding pairing interaction is attractive, and the instability to the d-wave superconductivity is favored. It is noted that the effective susceptibility of the s-wave channel increases rapidly and becomes positive near half filling. Both the s- and d-wave channels have finite values at half filling (ρ=1). It is known that superconductivity should not appear in the undoped system. The finite value can be understood as a possible tendency to the corresponding paring state. When the system is doped, the two channels behave differently. The d-wave is strengthened, while the s-wave is weakened. This implies that an instability to d-wave channel develops upon doping. We do not show the results for the filling , where the values of the d-wave are almost zero, and it is less possible to have a superconducting instability.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (a) The effective susceptibility of several pairing channels as a function of filling on L = 12 lattice. (b) Size dependence of the effective susceptibility of d-wave symmetry. The inset shows the corresponding average signs for the linear sizes L=8,10,12,16. In both figures, the temperature is T=1/12, and the Hubbard interaction is U=2t.

We have shown that the d-wave pairing is the leading superconducting instability on a finite-size system. It is necessary to check its size dependence. Figure 2(b) plots the effective susceptibility of the d-wave state for sizes L = 8, 10, 12, 16. While the size dependence is negligible on either side of the ρ axis, it is pronounced in the middle region, i.e., , surrounding the VHS. The sign problem is severe in this region, and the data have large error bars. Nevertheless, the average sign is still better than 0.7 at T=1/12 for U=2t. As we enlarge the lattice sizes, the sign problem is improved and the value of increases. The reason for this may be the finite-size effect on the density of states, which contains discrete peaks instead of being continuous on a finite lattice. For larger lattice sizes, there are more states near the VHS, which may result in an enhanced superconducting response.

It is also observed that the curve of the effective susceptibility is peaked at about at T=1/12. As the temperature is increased, the position of the peak shifts to larger ρ and finally locates at half filling (see Fig. 3(a)). This represents the process of weakening the superconducting state by increasing the temperature. It is expected that the peak moves leftwards when the temperature is lowered and the optimal filling of the ground state should be less than ρ=0.9.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The effective susceptibility of d-wave channel as a function of filling at several temperatures for (a) U=2t, (b) U =3t, (c) U=4t; (d) as a function of filling for several values of U at a fixed temperature. Insets in panels (b) and (c) show the corresponding average signs. In all figures, the linear lattice size is L = 12.

The effective susceptibility should be divergent at the superconducting transition temperature. However, due to the sign problem, the simulations are limited to relatively high temperatures, which are above the transition point. We can only deduce the possible low-temperature behavior based on the trend of the high-temperature data. Figure 3(a) shows the effective susceptibility of the d-wave channel as a function of density ρ at several temperatures for U=2t. As the temperature is lowered, is increased, implying that the superconducting instability is enhanced. From the trend of the T-dependent curve at a fixed density, it is possible that is divergent at a finite low temperature. We also perform simulations for stronger couplings, where the sign problem becomes severe. We make the temperatures as low as possible while maintaining acceptable average signs. As shown in Figs. 3(b) and 3(c), the average sign almost approaches zero near the VHS at T=1/10 for U=3t and T=1/6 for U=4t, which results in a huge error bar. To generate reliable results, the simulations have to be limited to higher temperatures.

To study the evolution of the dominating pairing symmetry with the interaction, we plot as a function of filling at T=1/4 for both weak and strong couplings. As shown in Fig. 3(d), all curves have the maximum values at half filling, implying that T=1/4 is much higher than the superconducting transition point. Nevertheless, the values monotonously increase as the couplings are increased for high-filling levels, implying that the d-wave state is enhanced by strong couplings. In the low-filling region, the value of is negative, which suggests that the superconductivity is further suppressed by strong couplings. Meanwhile, the values of the p-wave channel are still negative and decrease with increasing interactions. No signature of a transition to a dominating p-wave symmetry is observed upon increasing the interactions.^[22]

We present systematic DQMC results of the dominating pairing symmetry in the doped honeycomb lattice. The simulations cover a full range of fillings, and both weak and strong couplings are considered. For weak couplings, the d-wave state is dominant. From the evolution of the effective susceptibility with the temperature, the optimal filling is estimated to be around the VHS filling. Although the simulations are limited to high temperatures for strong couplings, the values of the effective susceptibility at high-filling levels increase as the interaction is strengthened, which implies that the d-wave state is enhanced. Our investigation extends the existing DQMC simulations, which only focus on specific fillings for weak couplings. Moreover, our DQMC simulations for strong couplings and on larger sizes provide new insights for the understanding of the superconducting pairing symmetry in the doped honeycomb lattice.